BMI Percentile by SES and Individual Differences in Adolescents

Data for this study come from the subset of responses collected on the SAPA-project.org website between February 17, 2017 and July 22, 2019. The initial date is the day that the semi-random presentation of items to participants was changed to increase presentation of SPI-135 items, which are the basis for personality measurement in this study. This period also represents a new period of data collection on SAPA containing data that are not available in the public domain at the time of analysis. The end date of data collection was the first day following preregistration of analysis that the authors were able to analyze data.

Clean data

set.seed(052319)

# load packages
packages = c("tidyverse", "janitor", "psych", "devtools", 
             "PAutilities", "measurements", "here", "caret")
lapply(packages, library, character.only = TRUE)
rm(packages)

#read in data
load(here("../../SAPA data/original data/SAPAdata07feb2017thru22jul2019forSara2.rdata"))
sapa = SAPAdata07feb2017thru22jul2019x

source(here("scripts/personality_scales.R"))
keys = read.csv("data/superKey.csv", header = TRUE, row.names = 1)

# super key -- this contains the master key list for all of SAPA. every item ever administered and every scale you can score
# each row is a single item
# each column is a scale
# the value of a cell is 0 if that item is not part of that scale, 1 if that item positively loads on the scale, and -1 if the item negatively loads on the scale

Participants were included in the analysis if they were under the age of 18, from the United States, and had reported their biological sex at birth, height, and weight.

# remove participants who are 18 years or older and from the US
sapa = sapa %>%
  filter(age < 18) %>%
  filter(country == "USA") %>%
  filter(!is.na(sex)) %>%
  filter(!is.na(height)) %>%
  filter(!is.na(weight)) %>%
  filter(!is.na(p1edu) | !is.na(p2edu) |
           !is.na(p1occIncomeEst) | !is.na(p2occIncomeEst) |
           !is.na(p1occPrestige) | !is.na(p2occPrestige))

Parental education was transformed into a numeric variable: 1 (less than 12 years of education), 2 (high school graduate or GED), 3 (some college), 4 (currently in college/university), 5 (Associate’s degree), 6 (university degree), 7 (currently in graduate or professional school), and 8 (graduate or professional degree). All parental SES variables – education, estimated income and estimated prestige, were standardized to the sample and averaged to create a single index of parental SES.

# make sure occupational variables are numeric
sapa = sapa %>%
  mutate_at(vars(matches("^(p)\\d(occ)")), as.numeric)

#or years
sapa = sapa %>%
  mutate(p1edu = case_when(
    p2edu == "less12yrs" ~ "6", 
    p2edu == "HSgrad" ~ "12", 
    p2edu == "SomeCollege" ~ "14", 
    p2edu == "CurrentInUniv" ~ "14",   
    p2edu == "AssociateDegree" ~ "14", 
    p2edu == "CollegeDegree" ~ "16", 
    p2edu == "InGradOrProSchool" ~ "18", 
    p2edu == "GradOrProDegree" ~ "20")) 

sapa = sapa %>%
  mutate(p2edu = case_when(
    p2edu == "less12yrs" ~ "6", 
    p2edu == "HSgrad" ~ "12", 
    p2edu == "SomeCollege" ~ "14", 
    p2edu == "CurrentInUniv" ~ "14",   
    p2edu == "AssociateDegree" ~ "14", 
    p2edu == "CollegeDegree" ~ "16", 
    p2edu == "InGradOrProSchool" ~ "18", 
    p2edu == "GradOrProDegree" ~ "20")) 

sapa$p1edu = as.numeric(sapa$p1edu)
sapa$p2edu = as.numeric(sapa$p2edu)

#estimate SES composite

sapa = sapa %>%
  mutate(z.p1edu = scale(p1edu),
         z.p2edu = scale(p2edu),
         z.p1occIncomeEst = scale(p1occIncomeEst),
         z.p2occIncomeEst = scale(p2occIncomeEst),
         z.p1occPrestige = scale(p1occPrestige),
         z.p2occPrestige = scale(p2occPrestige)) 

sapa$ses = rowMeans(sapa[,grepl("^z\\.", names(sapa))], na.rm=T)

sapa = sapa %>%
  dplyr::select(-starts_with("z"))

Big Five traits were scored using a sum-score method, averaged across non-missing responses.

# select just the rows that correspond to variables in the current SAPA dataset
vars = names(sapa)
keys = keys[rownames(keys) %in% vars, ]

# select just the Big 5 scales that are scored using the SPI_135 form 
bfkeys = keys %>%
  select(contains("SPI_135")) %>%
  select(1:5) 

bfkeys = keys2list(as.matrix(bfkeys), sign = T)


# score the items (this contains item and scale statistics too!)
b5scored = scoreItems(keys = bfkeys, items = sapa)

# add scores to SAPA
b5scores = as.data.frame(b5scored$scores[,1:5])
names(b5scores) = gsub("135_27_5_", "", names(b5scores))
sapa = cbind(sapa, b5scores)

The narrower traits, the SPI-27, were scored using IRT scoring. Calibration parameters were taken from a different dataset and are available on request.

load(here("../../SAPA data/created/IRTinfoSPI27.rdata"))

# IRT score
dataSet <- subset(sapa, select = c(orderForItems))

SPIirtScores <- matrix(nrow=dim(dataSet)[1], ncol=27)

scaleNames = gsub("SPI27_", "", names(IRToutputSPI27))
spi_keys = keys %>%
  select(matches("SPI_135")) %>%
  select(-c(1:5)) %>%
  mutate(item = rownames(.)) %>%
  gather("scale", "key", -item) %>%
  filter(key != 0)

for (i in 1:length(IRToutputSPI27)) {
  data <- subset(dataSet, select = c(rownames(IRToutputSPI27[[i]]$irt$difficulty[[1]])))
  calibrations <- IRToutputSPI27[[i]]
  #check calibration direction
  loadings = calibrations$fa$loadings[,1]
  loadings = ifelse(loadings < 0, -1, 1)
  loadings = data.frame(item = names(loadings), loadings = loadings)
  keys_direction = spi_keys %>%
    filter(grepl(scaleNames[i], scale)) %>%
    full_join(loadings)
  same = sum(keys_direction$key == keys_direction$loadings)
  if(same == 0) data[,1:ncol(data)] = apply(data[,1:ncol(data)], 2, function(x) max(x, na.rm=T) + 1 - x)
  if (same > 0 & same < 5) print("Error in loadings")
  scored <- scoreIrt(calibrations, data, keys = NULL, cut = 0)
  trait_scores = scored$theta1
  trait_scores = (trait_scores - mean(trait_scores, na.rm = T))/sd(trait_scores, na.rm=T)
  Tscores = trait_scores*10 + 50
  SPIirtScores[,i] <- Tscores
}

SPIirtScores <- as.data.frame(SPIirtScores)
colnames(SPIirtScores) <- paste0("SPI_", scaleNames)

#add to sapa dataset
sapa = cbind(sapa, SPIirtScores)

Cognition was also scored using IRT scoring, with calibrations from a separate dataset.

load(here("../../SAPA data/created/IRTinfoSPI27.rdata"))

# IRT score
dataSet <- subset(sapa, select = c(orderForItems))

SPIirtScores <- matrix(nrow=dim(dataSet)[1], ncol=27)

scaleNames = gsub("SPI27_", "", names(IRToutputSPI27))
spi_keys = keys %>%
  select(matches("SPI_135")) %>%
  select(-c(1:5)) %>%
  mutate(item = rownames(.)) %>%
  gather("scale", "key", -item) %>%
  filter(key != 0)

for (i in 1:length(IRToutputSPI27)) {
  data <- subset(dataSet, select = c(rownames(IRToutputSPI27[[i]]$irt$difficulty[[1]])))
  calibrations <- IRToutputSPI27[[i]]
  #check calibration direction
  loadings = calibrations$fa$loadings[,1]
  loadings = ifelse(loadings < 0, -1, 1)
  loadings = data.frame(item = names(loadings), loadings = loadings)
  keys_direction = spi_keys %>%
    filter(grepl(scaleNames[i], scale)) %>%
    full_join(loadings)
  same = sum(keys_direction$key == keys_direction$loadings)
  if(same == 0) data[,1:ncol(data)] = apply(data[,1:ncol(data)], 2, function(x) max(x, na.rm=T) + 1 - x)
  if (same > 0 & same < 5) print("Error in loadings")
  scored <- scoreIrt(calibrations, data, keys = NULL, cut = 0)
  trait_scores = scored$theta1
  trait_scores = (trait_scores - mean(trait_scores, na.rm = T))/sd(trait_scores, na.rm=T)
  Tscores = trait_scores*10 + 50
  SPIirtScores[,i] <- Tscores
}

SPIirtScores <- as.data.frame(SPIirtScores)
colnames(SPIirtScores) <- paste0("SPI_", scaleNames)

#add to sapa dataset
sapa = cbind(sapa, SPIirtScores)

BMI percentile represents a participant’s percentile score on BMI relative to others of their assigned sex at birth. These were estimated from the PAutilities package, developed by WHO Multicentre Growth Reference Study (MGRS) information about the development of these reference standards can be found at https://www.cdc.gov/obesity/childhood/defining.html. These standards in turn were develoed using the 2000 CDC growth charts, based on data from 5 national health examination surveys that occurred from 1963 to 1994 and supplemental data from surveys that occurred from 1960 to 1995.

Kuczmarski RJ, Ogden CL, Guo SS, et al. 2000 CDC growth charts for the United States: methods and development. National Center for Health Statistics. Vital Health Stat 11. 2002;(246):1-190

BMI category is assigned based on BMI percentile: participants in the bottom 10% are labeled Underweight, between the top 10% and 5% are Overweight, and top 5% are Obese. All others are labelled Normal.

All analyses were perfomed separately by gender.

sapa = sapa %>%
  mutate(cog = ICAR60) %>%
  select(sex, age, height, weight, BMI, BMI_p, BMI_c, p1edu, 
         p1occPrestige, p1occIncomeEst, p2edu, 
         p2occPrestige, p2occIncomeEst, ses, cog, contains("SPI"))

sapa_male = sapa %>%
  filter(sex == "male") %>%
  dplyr::select(-sex) 

sapa_female = sapa %>%
  filter(sex == "female") %>%
  dplyr::select(-sex)

save(b5scored, file = here("data/alpha.Rdata"))

The datasets were split into training (75%) and test (25%) sets; all regression models are estimated using the training sets. The test set was reserved to estimate model accuarcy, comparing models with different sets of individual differences.

# set seed
set.seed(090919)

# parition into training and test sets. objects identify just training rows
train_male = createDataPartition(sapa_male$BMI_c, p = .75, list = FALSE)
train_female = createDataPartition(sapa_female$BMI_c, p = .75, list = FALSE)

Descriptive Statistics

Univariate descriptives

Descriptive statistics are estimated using the psych package.

descriptives = describeBy(sapa, group = "sex")

library(kableExtra)
#pull descriptives into a list
descriptives.df = data.frame(gender = names(descriptives))
descriptives.df$data = descriptives

#add variable names and unnest
descriptives.df = descriptives.df %>%
  mutate(data = map(data, function(x) mutate(x, vars = rownames(x)))) %>%
  unnest(cols = c(data))

Descriptives Table by Gender

Female

descriptives.df %>%
  filter(gender == "female") %>%
  select(-gender) %>%
  kable(., digits = 2) %>%
  kable_styling()

vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
sex*	6530	1.00	0.00	1.00	1.00	0.00	1.00	1.00	0.00	NaN	NaN	0.00
BMI	6502	23.07	5.00	21.85	22.38	3.65	15.01	52.78	37.77	1.65	3.85	0.06
BMI_p	6530	62.70	27.61	67.80	64.88	30.99	0.00	99.80	99.80	-0.53	-0.80	0.34
BMI_c*	6530	1.51	0.84	1.00	1.35	0.00	1.00	4.00	3.00	1.38	0.54	0.01
p1edu	5876	14.62	3.66	14.00	14.87	2.97	6.00	20.00	14.00	-0.53	0.36	0.05
p1occPrestige	5723	60.76	14.64	67.12	62.14	13.55	24.22	79.09	54.87	-0.82	-0.42	0.19
p1occIncomeEst	5594	61625.23	21784.89	60244.00	60512.80	20228.59	21980.00	112490.00	90510.00	0.39	-0.24	291.27
p2edu	5876	14.62	3.66	14.00	14.87	2.97	6.00	20.00	14.00	-0.53	0.36	0.05
p2occPrestige	4818	57.87	15.76	63.85	58.75	16.23	24.22	79.09	54.87	-0.47	-1.04	0.23
p2occIncomeEst	4729	59058.07	22926.91	55790.00	57519.07	24640.81	21980.00	112490.00	90510.00	0.53	-0.34	333.40
ses	6452	-0.03	0.78	0.11	0.02	0.75	-2.38	1.55	3.93	-0.63	0.22	0.01
cog	6507	48.06	9.19	48.02	48.00	9.73	18.46	76.53	58.07	0.03	-0.35	0.11
SPI_Agree	6530	4.26	0.67	4.36	4.31	0.53	1.00	6.00	5.00	-0.81	1.35	0.01
SPI_Consc	6530	3.87	0.67	3.86	3.87	0.64	1.00	6.00	5.00	-0.05	0.36	0.01
SPI_Extra	6530	3.61	0.83	3.64	3.62	0.85	1.00	6.00	5.00	-0.09	-0.02	0.01
SPI_Neuro	6530	4.31	0.77	4.36	4.33	0.74	1.29	6.00	4.71	-0.42	0.28	0.01
SPI_Open	6530	4.48	0.55	4.50	4.49	0.53	2.07	6.00	3.93	-0.24	0.31	0.01
SPI_Compassion	6530	51.31	9.51	53.78	52.56	8.67	15.93	63.13	47.20	-1.15	1.12	0.12
SPI_Irritability	6530	50.97	9.87	51.27	51.15	11.81	27.87	72.50	44.63	-0.13	-0.96	0.12
SPI_Sociability	6530	49.82	10.02	52.22	51.19	8.27	15.89	66.44	50.54	-1.36	1.94	0.12
SPI_WellBeing	6488	49.03	9.80	48.95	49.04	11.22	26.47	74.40	47.94	-0.01	-0.84	0.12
SPI_SensationSeeking	6488	49.55	10.03	49.00	49.26	11.45	29.19	71.62	42.43	0.21	-0.88	0.12
SPI_Anxiety	6530	51.88	8.90	54.46	53.24	7.80	19.27	62.10	42.82	-1.26	1.19	0.11
SPI_Honesty	6530	50.65	9.56	52.54	51.80	7.38	3.87	77.00	73.13	-1.70	4.87	0.12
SPI_Industry	6530	50.29	10.01	49.64	50.18	11.27	27.91	73.23	45.31	0.10	-0.80	0.12
SPI_Intellect	6530	49.36	10.20	51.63	50.56	9.26	14.10	65.73	51.63	-1.04	0.75	0.13
SPI_Creativity	6530	49.72	10.03	51.03	50.57	10.36	18.38	64.67	46.29	-0.70	0.03	0.12
SPI_Impulsivity	6488	49.85	10.09	49.16	49.47	11.50	31.84	72.72	40.88	0.27	-0.82	0.13
SPI_AttentionSeeking	6530	49.68	10.08	52.47	50.98	8.25	26.43	66.95	40.51	-1.00	0.12	0.12
SPI_Order	6529	50.18	10.06	49.75	49.98	11.63	25.20	74.44	49.24	0.15	-0.87	0.12
SPI_Authoritarianism	6489	50.55	9.69	52.46	51.61	8.97	15.22	67.09	51.87	-0.95	0.52	0.12
SPI_Charisma	6487	49.74	10.00	50.32	50.11	10.90	22.58	71.58	49.00	-0.30	-0.57	0.12
SPI_Trust	6530	49.83	10.00	50.04	49.87	10.68	27.80	73.20	45.40	-0.03	-0.71	0.12
SPI_Humor	6489	50.65	9.62	52.72	51.87	8.77	7.52	64.74	57.21	-1.25	1.81	0.12
SPI_EmotionalExpressiveness	6489	50.36	10.04	50.04	50.12	12.09	32.70	71.81	39.11	0.16	-1.01	0.12
SPI_ArtAppreciation	6488	51.35	8.94	53.85	52.92	5.09	15.08	72.72	57.65	-2.11	5.33	0.11
SPI_Introspection	6529	50.04	9.89	51.98	51.19	9.80	15.16	62.58	47.43	-0.95	0.48	0.12
SPI_Perfectionism	6530	50.73	9.83	51.36	51.18	10.47	19.29	69.55	50.25	-0.38	-0.44	0.12
SPI_SelfControl	6485	49.31	9.86	49.50	49.30	10.53	25.93	74.34	48.40	0.01	-0.56	0.12
SPI_Conformity	6489	50.59	9.93	51.34	50.83	10.79	28.22	71.14	42.92	-0.20	-0.77	0.12
SPI_Adaptability	6488	49.44	10.12	48.75	49.39	11.32	28.77	70.29	41.52	0.06	-0.88	0.13
SPI_EasyGoingness	6485	49.56	9.99	50.80	50.22	9.63	13.49	69.53	56.03	-0.61	0.09	0.12
SPI_EmotionalStability	6529	48.20	10.02	50.26	48.89	9.90	28.98	68.44	39.46	-0.56	-0.68	0.12
SPI_Conservatism	6480	49.54	10.13	51.27	50.42	10.20	28.73	70.97	42.24	-0.63	-0.46	0.13

Male

descriptives.df %>%
  filter(gender == "male") %>%
  select(-gender) %>%
  kable(., digits = 2) %>%
  kable_styling()

vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
sex*	2952	2.00	0.00	2.00	2.00	0.00	2.00	2.00	0.00	NaN	NaN	0.00
BMI	2934	22.84	4.90	21.62	22.21	3.69	15.01	53.76	38.75	1.48	3.05	0.09
BMI_p	2952	60.00	30.53	64.60	62.02	37.95	0.00	99.90	99.90	-0.40	-1.07	0.56
BMI_c*	2952	1.57	0.91	1.00	1.40	0.00	1.00	4.00	3.00	1.37	0.61	0.02
p1edu	2660	14.71	3.58	16.00	14.95	2.97	6.00	20.00	14.00	-0.60	0.55	0.07
p1occPrestige	2620	60.20	15.22	67.12	61.55	13.55	24.22	79.09	54.87	-0.75	-0.62	0.30
p1occIncomeEst	2570	61491.45	22195.84	60244.00	60663.71	22141.15	21980.00	112490.00	90510.00	0.29	-0.44	437.83
p2edu	2660	14.71	3.58	16.00	14.95	2.97	6.00	20.00	14.00	-0.60	0.55	0.07
p2occPrestige	2147	57.07	15.59	57.98	58.02	18.55	24.22	79.09	54.87	-0.50	-0.98	0.34
p2occIncomeEst	2111	57247.11	22364.35	52210.00	55755.02	24789.07	21980.00	112490.00	90510.00	0.60	-0.04	486.76
ses	2917	-0.05	0.79	0.05	0.00	0.77	-2.38	1.55	3.93	-0.64	0.19	0.01
cog	2948	49.65	9.60	49.71	49.68	10.12	19.86	76.53	56.68	-0.06	-0.37	0.18
SPI_Agree	2952	4.11	0.70	4.21	4.17	0.53	1.00	6.00	5.00	-1.00	1.88	0.01
SPI_Consc	2952	3.77	0.63	3.79	3.77	0.53	1.43	6.00	4.57	-0.01	0.64	0.01
SPI_Extra	2952	3.62	0.81	3.64	3.63	0.74	1.00	6.00	5.00	-0.20	0.24	0.01
SPI_Neuro	2952	3.79	0.81	3.86	3.81	0.74	1.00	6.00	5.00	-0.27	0.29	0.01
SPI_Open	2952	4.55	0.53	4.57	4.56	0.53	2.21	6.00	3.79	-0.27	0.50	0.01
SPI_Compassion	2952	47.10	10.40	48.58	48.07	9.76	15.93	63.13	47.20	-0.83	0.40	0.19
SPI_Irritability	2952	47.83	9.94	47.40	47.58	11.58	27.87	72.50	44.63	0.18	-0.91	0.18
SPI_Sociability	2952	50.46	9.87	53.03	51.88	7.51	15.89	66.44	50.54	-1.51	2.53	0.18
SPI_WellBeing	2925	52.24	10.04	53.13	52.67	11.35	26.47	74.40	47.94	-0.30	-0.77	0.19
SPI_SensationSeeking	2921	50.98	9.86	50.67	50.85	10.93	30.74	71.42	40.68	0.11	-0.88	0.18
SPI_Anxiety	2952	45.79	10.97	47.72	46.64	11.56	17.73	62.10	44.37	-0.58	-0.59	0.20
SPI_Honesty	2952	48.62	10.74	50.53	49.85	8.06	3.87	77.00	73.13	-1.49	3.79	0.20
SPI_Industry	2952	49.39	9.93	48.86	49.13	10.98	28.29	73.92	45.64	0.20	-0.73	0.18
SPI_Intellect	2952	51.43	9.39	53.71	52.62	8.79	14.10	65.73	51.63	-1.15	1.25	0.17
SPI_Creativity	2952	50.56	9.92	52.20	51.47	10.46	18.38	64.67	46.29	-0.77	0.17	0.18
SPI_Impulsivity	2921	50.31	9.78	49.82	49.99	10.76	31.84	72.72	40.88	0.24	-0.76	0.18
SPI_AttentionSeeking	2952	50.72	9.77	53.52	52.13	7.58	26.43	65.70	39.27	-1.14	0.54	0.18
SPI_Order	2951	49.64	9.87	49.08	49.39	11.12	25.20	74.44	49.24	0.18	-0.77	0.18
SPI_Authoritarianism	2925	48.87	10.55	51.21	49.93	9.93	13.46	67.09	53.63	-0.86	0.23	0.20
SPI_Charisma	2926	50.59	9.97	51.49	51.01	10.96	21.90	71.58	49.68	-0.34	-0.54	0.18
SPI_Trust	2952	50.44	9.94	50.96	50.57	10.52	27.80	73.20	45.40	-0.10	-0.61	0.18
SPI_Humor	2923	48.57	10.62	50.33	49.73	10.16	7.52	64.74	57.21	-1.02	1.04	0.20
SPI_EmotionalExpressiveness	2923	49.18	9.84	48.23	48.76	11.20	32.70	71.81	39.11	0.31	-0.88	0.18
SPI_ArtAppreciation	2921	46.96	11.47	50.20	48.80	7.78	15.08	72.72	57.65	-1.35	1.37	0.21
SPI_Introspection	2951	49.92	10.21	52.19	51.08	10.56	15.16	62.58	47.43	-0.88	0.19	0.19
SPI_Perfectionism	2952	48.37	10.18	48.80	48.70	10.96	18.72	69.55	50.82	-0.27	-0.47	0.19
SPI_SelfControl	2919	51.54	10.12	51.77	51.70	10.72	25.93	74.34	48.40	-0.13	-0.51	0.19
SPI_Conformity	2926	48.76	10.03	49.07	48.80	11.21	28.22	70.64	42.42	-0.03	-0.82	0.19
SPI_Adaptability	2921	51.24	9.60	51.22	51.44	10.68	28.77	70.29	41.52	-0.13	-0.75	0.18
SPI_EasyGoingness	2920	50.97	9.98	52.52	51.78	9.37	13.49	69.53	56.03	-0.75	0.33	0.18
SPI_EmotionalStability	2952	54.01	8.69	56.24	55.20	6.98	28.98	67.27	38.29	-1.22	1.17	0.16
SPI_Conservatism	2919	51.07	9.62	52.89	52.08	9.05	28.73	70.97	42.24	-0.79	-0.02	0.18

Bivariate

R_male = sapa_male %>%
  dplyr::select(-BMI_c) %>%
  cor(use = "pairwise")

R_female = sapa_female %>%
  dplyr::select(-BMI_c) %>%
  cor(use = "pairwise")



#predictors
pred = names(sapa_male) %>% str_subset("BMI", negate = TRUE)


r_bmi_male = corr.test(x = sapa_male$BMI, y = sapa_male[,pred])
r_bmi_female = corr.test(x = sapa_female$BMI, y = sapa_female[,pred])

r_bmi_male = modify(r_bmi_male, as.vector)
r_bmi_female = modify(r_bmi_female, as.vector)

cor.data = data.frame(gender = c("male", "female"))
cor.data$fullr = list(r_bmi_male, r_bmi_female) 


cor.data = cor.data %>%
  mutate(r = map(fullr, "r")) %>%
  mutate(r = map(r, unlist)) %>%
  mutate(rp = map(fullr, "p")) %>%
  mutate(rp = map(rp, unlist)) %>%
  dplyr::select(-fullr) %>%
  unnest(cols = c(r, rp)) %>%
  mutate(pred = rep(pred,2)) %>%
  gather("key", "value", -gender, -pred) %>%
  unite(gender, gender, key) %>%
  spread(gender, value)

save(R_male, R_female, cor.data, file = "data/cor_output.Rdata")

Figure

Female

library(corrplot)
corrplot(R_female, method = "square",
         title = "\nZero-order correlations among study variables\nFemale Participants",
         tl.col = "black",
         mar=c(0,0,1,0))

Male

corrplot(R_male, method = "square",
         title = "\nZero-order correlations among study variables\nMale Participants",
         tl.col = "black",
         mar=c(0,0,1,0))

Regression Models (SES)

Regression models were built that regressed BMI percentile onto parental socio-economic status and adolescent individual differences. Two basic models were constructed: one that hypothesized parental SES:

\[BMIp_i = b_0 + b_1(SES_i) + b_2(ID_i) + e_i\]

and an individual difference were two independent predictors of BMI, and a second that hypothesized these variables interacted with each other:

\[BMIp_i = b_0 + b_1(SES_i) + b_2(ID_i) + b_3(SES_i\times ID_i) + e_i\] We iterated through all individual differences – the broad Big Five personality traits, the narrow SPI-27 traits, and cognitive ability – and tested each one independently in the model as an individual difference.

Models were estimated separately for men and women.

#end goal of wrangling is a data frame of data frames
# nested dataframes correspond to a single personality trait
# score refers to a participant's score on that trait
# we also standardize each of our variables within gender


sapa_male_trait = sapa_male[train_male, ] %>%
  dplyr::select(-starts_with("p1"), -starts_with("p2")) %>%
  mutate(set = ifelse(row_number() %in% train_male[,1], "train", "test")) %>% #identify which rows in test and training
  gather("trait_name", "trait_score", -ses, -BMI_c, -BMI, -BMI_p, -set) %>% # gather all personality variables
  group_by(trait_name, set) %>% # group by trait and also by whether in test/train
  mutate(trait_score = scale(trait_score)) %>% #standardize
  mutate(ses = scale(ses)) %>% #standardize
  ungroup() %>% group_by(trait_name) %>%  #group only by trait
  nest() #nest data frames


sapa_female_trait = sapa_female[train_female, ] %>%
  dplyr::select(-starts_with("p1"), -starts_with("p2")) %>%
  mutate(set = ifelse(row_number() %in% train_male[,1], "train", "test")) %>% #identify which rows in test and training
  gather("trait_name", "trait_score", -ses, -BMI_c, -BMI, -BMI_p, -set) %>% # gather all personality variables
  group_by(trait_name, set) %>% # group by trait and also by whether in test/train
  mutate(trait_score = scale(trait_score)) %>% #standardize
  mutate(ses = scale(ses)) %>% #standardize
  ungroup() %>% group_by(trait_name) %>%  #group only by trait
  nest() #nest data frames

Controlling for personality

To estimate the effect of socioeconomic status on BMI percentile, we graph the estimates of the SES slope coefficient across all regression models controlling for individual differences. This presents not only the average estimate across all models (solid line), but the range of estimates – a wide range suggests that the effect of SES on BMI is sensitive to the inclusion of different individual difference measures, while a narrow range suggests that the effect of SES on BMI is persistent through personality and cognition.

Female

avg_female = female_reg %>%
  filter(term == "ses") %>%
  filter(model == "cov") %>%
  summarize(mean = mean(estimate))

female_plot_1 = female_reg %>%
  filter(term == "ses") %>%
  filter(model == "cov") %>%
  mutate(psig = ifelse(p.value < .05, "yes", "no")) %>%
  arrange(estimate) %>%
  mutate(spec = row_number()) %>%
  ggplot(aes(x = spec, y = conf.low)) +
  geom_segment(aes(xend = spec, yend = conf.high, color = psig)) + 
  #geom_point(aes(y = estimate, color = "grey")) +
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  geom_hline(aes(yintercept = mean), data = avg_female) +
  #geom_label(aes(x = 25, y = 1.25, label = round(mean,2)), data = avg_female )+
  scale_color_manual(values = c("red", "grey")) +
  scale_y_continuous(limits = c(-5.5, 0.25), breaks = c(-5:0))+
  labs(x = "Specification",
       y = "SES coefficient, controlling for personality", 
       title = "Adolescent Girls") +
  guides(color = F) +
  theme_pubr()

female_plot_1

Male

avg_male = male_reg %>%
  filter(term == "ses") %>%
  filter(model == "cov") %>%
  summarize(mean = mean(estimate))

male_plot_1 = male_reg %>%
  filter(term == "ses") %>%
  filter(model == "cov") %>%
  mutate(psig = ifelse(p.value < .05, "yes", "no")) %>%
  arrange(estimate) %>%
  mutate(spec = row_number()) %>%
  ggplot(aes(x = spec, y = conf.low)) +
  geom_segment(aes(xend = spec, yend = conf.high, color = psig)) + 
  #geom_point(aes(y = estimate)) +
  geom_hline(aes(yintercept = 0), linetype = "dashed", color = "black") +
  geom_hline(aes(yintercept = mean), data = avg_male, color = "black") +
  #geom_label(aes(x = 25, y = 1.25, label = round(mean,2)), data = avg_male )+
  scale_color_manual(values = c("red", "grey")) +
  scale_y_continuous(limits = c(-5.5, .25), breaks = c(-5:0))+
  labs(x = "Specification", title = "Adolescent Boys", y = NULL) +
  guides(color = F) +
  theme_pubr()

male_plot_1

Interaction with personality

To estimate the joint effect of socioeconomic status and individual differences on BMI percentile, we graph the estimates of the interaction terms of SES by individual differences by BMI percentile. Like before, we present the average effect (solid black line) and the 95% confidence intervals for each model.

Female

avg_female = female_reg %>%
  filter(grepl(":", term)) %>%
  summarize(mean = mean(estimate))

Parental socioeconomic status positively predicted greater likelihood of all non-normal categories (Underweight, Overweight, and Obese) compared to Normal among girls. Adolescent girls living in higher SES households were, on average, NA% less likely to be Underweight, NA% less likely, and 97% less likely to be Obese compared to low SES counterparts.

female_plot_2 = female_reg %>%
  filter(grepl(":", term)) %>%
  mutate(psig = ifelse(p.value < .05, "yes", "no")) %>%
  arrange(estimate) %>%
  mutate(spec = row_number()) %>%
  ggplot(aes(x = spec, y = conf.low)) +
  geom_segment(aes(xend = spec, yend = conf.high, color = psig)) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed", color = "black") +
  geom_hline(aes(yintercept = mean), data = avg_female, color = "black") +
  #geom_label(aes(x = 25, y = 1.25, label = round(mean,2)), data = avg_female )+
  scale_color_manual(values = c("grey", "red")) +
  #scale_y_continuous(limits = c(0.30, 2.20))+
  labs(x = "Specification",
       y = "SES x perosnality term in model", title = "Adolescent Girls") +
  guides(color = F) +
  theme_pubr()

female_plot_2

Male

avg_male = male_reg %>%
  filter(grepl(":", term)) %>%
  summarize(mean = mean(estimate))

Parental socioeconomic status positively predicted greater likelihood of all non-normal categories (Underweight, Overweight, and Obese) compared to Normal among boys. Adolescent boys living in higher SES households were, on average, NA% more likely to be Underweight, NA% more likely, and -59% more likely to be Obese compared to low SES counterparts.

male_plot_2 = male_reg %>%
  filter(grepl(":", term)) %>%
  mutate(psig = ifelse(p.value < .05, "yes", "no")) %>%
  arrange(estimate) %>%
  mutate(spec = row_number()) %>%
  ggplot(aes(x = spec, y = conf.low)) +
  geom_segment(aes(xend = spec, yend = conf.high, color = psig)) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed", color = "black") +
  geom_hline(aes(yintercept = mean), data = avg_male, color = "black") +
  #geom_label(aes(x = 25, y = 1.25, label = round(mean,2)), data = avg_male )+
  scale_color_manual(values = c("grey", "red")) +
  #scale_y_continuous(limits = c(-5.5, .25))+
  labs(x = "Specification",
       y = "Odds Ratio of SES x perosnality term in model") +
  guides(color = F) +
  theme_pubr()

male_plot_2

# both plots together -----------------------------------------------------

ggarrange(female_plot_1, male_plot_1)

ggsave(here("figures/SES_specification.jpeg"), width = 6, height = 4)

Code (Female)

# we run the models for men and women separately because R kept crashing whey trying to run this whole script.2

female_reg = sapa_female_trait %>%
  mutate(cov = map(data, ~lm(BMI_p ~ trait_score + ses, data = .))) %>%
  mutate(int = map(data, ~lm(BMI_p ~ trait_score*ses, data = .))) 

female_plot = female_reg %>%
  mutate(cov = map(data, ~lm(BMI_p ~ trait_score + ses, data = .))) %>%
  mutate(int = map(data, ~lm(BMI_p ~ trait_score*ses, data = .))) 


female_reg = female_reg %>%
  dplyr::select(-data) %>%
  gather("model", "output", cov, int) %>%
  mutate(output = map(output, broom::tidy, conf.int = FALSE)) %>%
  unnest(cols = c(output))

Code (Male)

male_reg = sapa_male_trait %>%
  mutate(cov = map(data, ~lm(BMI_p ~ trait_score + ses, data = .))) %>%
  mutate(int = map(data, ~lm(BMI_p ~ trait_score*ses, data = .))) 

male_plot = male_reg %>%
  mutate(cov = map(data, ~lm(BMI_p ~ trait_score + ses, data = .))) %>%
  mutate(int = map(data, ~lm(BMI_p ~ trait_score*ses, data = .))) 


male_reg = male_reg %>%
  dplyr::select(-data) %>%
  gather("model", "output", cov, int) %>%
  mutate(output = map(output, broom::tidy, conf.int = FALSE)) %>%
  unnest(cols = c(output))

Regression (Individual differences)

load(here("data/regression_output_male.Rdata"))
load(here("data/regression_output_female.Rdata"))


names(SPI_5_names) = str_remove(names(SPI_5_names), "135_27_5_")
names(SPI_27_names) = str_remove(names(SPI_27_names), "135_27_5_")

female_reg = female_reg %>%
  mutate(gender = "female") %>%
  ungroup()

male_reg = male_reg %>%
  mutate(gender = "male") %>%
  ungroup()

female_reg %>%
  full_join(male_reg) %>%
  filter(grepl("trait", term)) %>%
  mutate(b1_est = printnum(estimate),  
         b1_est = ifelse(conf.low > 0 | conf.high < 0, paste0(b1_est, "*"), b1_est),
         conf.low = printnum(conf.low),
         conf.high = printnum(conf.high), 
         b2_conf = paste0("[", conf.low, ", ", conf.high, "]")) %>%
  dplyr::select(trait_name, model, term, b1_est, b2_conf, gender) %>%
  gather("key", "value", b1_est, b2_conf) %>%
  unite(col = "newkey", gender, model, term) %>%
  spread(newkey, value) %>%
  mutate(trait_name = factor(trait_name, 
                             levels = c("cog", names(SPI_27_names), names(SPI_5_names)),
                             labels = c("Cognitive Ability", SPI_27_names, SPI_5_names))) %>%
  arrange(trait_name) %>% 
  #mutate(trait_name = ifelse(!(row_number() %% 2), NA_character_, trait_name)) %>%
  dplyr::select(-key) %>%
  kable(., col.names = rep(c("Trait", rep(c("b", "b", "b x SES"), 2)))) %>%
  kable_styling() %>%
  add_header_above(c(" ", "Additive Model" = 1, "Interaction Model" = 2, "Additive Model" = 1, "Interaction Model" = 2)) %>%
  add_header_above(c(" ", "Female" = 3, "Male" = 3)) %>%
  group_rows("SPI: 27 Factors", 3, 56) %>%
  group_rows("SPI: 5 Factors", 57, 66) 

## Joining, by = c("trait_name", "model", "term", "estimate", "std.error", "statistic", "p.value", "conf.low", "conf.high", "gender")

## Warning: attributes are not identical across measure variables;
## they will be dropped

	Female			Male
	Additive Model	Interaction Model		Additive Model	Interaction Model
Trait	b	b	b x SES	b	b	b x SES
Cognitive Ability	-1.32*	-1.33*	0.09	-1.84*	-1.83*	-0.08
Cognitive Ability	[-2.10, -0.54]	[-2.11, -0.54]	[-0.68, 0.86]	[-3.10, -0.53]	[-3.10, -0.52]	[-1.37, 1.20]
SPI: 27 Factors
Compassion	-0.20	-0.19	-0.38	-0.31	-0.34	0.44
Compassion	[-0.98, 0.61]	[-0.98, 0.61]	[-1.14, 0.37]	[-1.61, 1.01]	[-1.64, 0.98]	[-0.81, 1.65]
Irritability	1.42*	1.43*	0.24	1.03	1.03	0.29
Irritability	[0.64, 2.20]	[0.65, 2.20]	[-0.49, 0.97]	[-0.26, 2.30]	[-0.26, 2.31]	[-0.98, 1.60]
Sociability	-1.21*	-1.21*	0.33	0.31	0.39	1.22
Sociability	[-2.01, -0.41]	[-2.01, -0.41]	[-0.47, 1.11]	[-0.97, 1.58]	[-0.88, 1.69]	[-0.04, 2.46]
Well Being	-2.72*	-2.70*	0.81*	-0.20	-0.18	0.57
Well Being	[-3.50, -1.94]	[-3.48, -1.93]	[0.04, 1.58]	[-1.53, 1.14]	[-1.52, 1.15]	[-0.69, 1.78]
Sensation Seeking	-0.58	-0.56	0.63	-0.03	-0.03	-0.40
Sensation Seeking	[-1.38, 0.22]	[-1.37, 0.24]	[-0.13, 1.41]	[-1.33, 1.25]	[-1.35, 1.25]	[-1.68, 0.90]
Anxiety	1.34*	1.38*	-0.50	-0.34	-0.35	0.65
Anxiety	[0.54, 2.14]	[0.59, 2.18]	[-1.29, 0.29]	[-1.61, 0.93]	[-1.62, 0.92]	[-0.60, 1.92]
Honesty	-1.03*	-1.04*	0.49	-0.19	-0.24	0.81
Honesty	[-1.78, -0.28]	[-1.79, -0.29]	[-0.24, 1.25]	[-1.43, 1.10]	[-1.47, 1.07]	[-0.40, 2.00]
Industry	-0.81*	-0.81*	-0.21	0.77	0.75	0.35
Industry	[-1.61, -0.04]	[-1.61, -0.04]	[-0.97, 0.54]	[-0.49, 2.05]	[-0.50, 2.03]	[-0.96, 1.61]
Intellect	-0.44	-0.45	-0.22	0.27	0.22	-0.55
Intellect	[-1.24, 0.34]	[-1.26, 0.33]	[-0.95, 0.51]	[-1.05, 1.55]	[-1.10, 1.49]	[-1.87, 0.80]
Creativity	-0.27	-0.27	0.02	0.22	0.22	0.11
Creativity	[-1.06, 0.51]	[-1.06, 0.51]	[-0.76, 0.77]	[-1.10, 1.52]	[-1.10, 1.52]	[-1.28, 1.53]
Impulsivity	0.78	0.77	0.39	0.00	0.01	-0.65
Impulsivity	[-0.03, 1.57]	[-0.04, 1.56]	[-0.42, 1.20]	[-1.28, 1.30]	[-1.26, 1.32]	[-1.98, 0.65]
Attention Seeking	-0.69	-0.65	0.50	-0.12	0.01	1.26
Attention Seeking	[-1.47, 0.11]	[-1.44, 0.15]	[-0.25, 1.25]	[-1.46, 1.21]	[-1.32, 1.35]	[-0.04, 2.55]
Order	-2.27*	-2.26*	-0.80*	-0.61	-0.60	-0.50
Order	[-3.03, -1.52]	[-3.02, -1.51]	[-1.54, -0.06]	[-1.92, 0.67]	[-1.90, 0.69]	[-1.81, 0.78]
Authoritarianism	0.37	0.37	0.17	0.52	0.44	1.51*
Authoritarianism	[-0.43, 1.17]	[-0.44, 1.17]	[-0.61, 0.96]	[-0.72, 1.78]	[-0.81, 1.68]	[0.25, 2.76]
Charisma	0.41	0.41	0.19	1.04	1.04	0.49
Charisma	[-0.38, 1.20]	[-0.38, 1.20]	[-0.56, 0.94]	[-0.24, 2.39]	[-0.24, 2.38]	[-0.81, 1.75]
Trust	-0.28	-0.28	0.02	-0.31	-0.40	0.96
Trust	[-1.06, 0.51]	[-1.06, 0.50]	[-0.77, 0.80]	[-1.60, 0.97]	[-1.68, 0.90]	[-0.29, 2.21]
Humor	1.03*	1.03*	-0.30	0.66	0.66	0.66
Humor	[0.23, 1.84]	[0.22, 1.84]	[-1.04, 0.44]	[-0.63, 1.96]	[-0.63, 1.96]	[-0.70, 2.02]
Emotional Expressiveness	-0.62	-0.63	0.33	-0.46	-0.53	1.36*
Emotional Expressiveness	[-1.41, 0.16]	[-1.42, 0.16]	[-0.46, 1.09]	[-1.78, 0.78]	[-1.84, 0.73]	[0.06, 2.66]
Art Appreciation	0.00	0.00	-0.19	-0.33	-0.33	-0.05
Art Appreciation	[-0.75, 0.73]	[-0.75, 0.74]	[-0.95, 0.55]	[-1.60, 0.94]	[-1.60, 0.94]	[-1.36, 1.19]
Introspection	-1.05*	-1.05*	0.37	-0.39	-0.37	0.47
Introspection	[-1.80, -0.28]	[-1.81, -0.29]	[-0.37, 1.08]	[-1.69, 0.91]	[-1.66, 0.92]	[-0.74, 1.69]
Perfectionism	-0.60	-0.61	-0.58	-0.93	-0.93	0.60
Perfectionism	[-1.40, 0.19]	[-1.41, 0.17]	[-1.33, 0.20]	[-2.18, 0.33]	[-2.18, 0.33]	[-0.66, 1.83]
Self Control	-2.79*	-2.79*	-0.07	-1.94*	-1.98*	1.00
Self Control	[-3.57, -1.99]	[-3.57, -1.99]	[-0.81, 0.67]	[-3.22, -0.65]	[-3.26, -0.70]	[-0.31, 2.34]
Conformity	0.90*	0.89*	-0.24	0.46	0.45	-0.19
Conformity	[0.10, 1.70]	[0.09, 1.70]	[-1.01, 0.55]	[-0.85, 1.77]	[-0.86, 1.76]	[-1.48, 1.03]
Adaptability	0.19	0.19	0.23	-0.40	-0.44	0.96
Adaptability	[-0.58, 0.94]	[-0.58, 0.94]	[-0.50, 0.97]	[-1.72, 0.90]	[-1.76, 0.87]	[-0.36, 2.29]
Easy Goingness	1.57*	1.59*	-0.33	1.09	1.19	-1.41*
Easy Goingness	[0.81, 2.35]	[0.82, 2.37]	[-1.11, 0.41]	[-0.17, 2.34]	[-0.08, 2.45]	[-2.67, -0.18]
Emotional Stability	-0.35	-0.35	0.23	1.00	1.00	-0.49
Emotional Stability	[-1.13, 0.45]	[-1.14, 0.45]	[-0.55, 1.01]	[-0.33, 2.29]	[-0.33, 2.29]	[-1.73, 0.73]
Conservatism	-0.94*	-0.97*	0.86*	1.32*	1.25	1.44*
Conservatism	[-1.72, -0.16]	[-1.77, -0.19]	[0.05, 1.65]	[0.01, 2.65]	[-0.05, 2.58]	[0.10, 2.83]
SPI: 5 Factors
Agreeableness	-0.36	-0.36	-0.28	-0.28	-0.37	0.76
Agreeableness	[-1.14, 0.41]	[-1.13, 0.41]	[-1.06, 0.52]	[-1.56, 1.02]	[-1.65, 0.94]	[-0.50, 2.01]
Conscientiousness	-1.35*	-1.33*	-0.76	-0.28	-0.28	0.49
Conscientiousness	[-2.12, -0.57]	[-2.10, -0.55]	[-1.54, 0.05]	[-1.58, 1.02]	[-1.58, 1.02]	[-0.73, 1.67]
Extraversion	-1.04*	-1.06*	0.56	0.51	0.56	1.45*
Extraversion	[-1.85, -0.24]	[-1.87, -0.27]	[-0.20, 1.31]	[-0.80, 1.80]	[-0.76, 1.86]	[0.13, 2.72]
Neuroticism	1.75*	1.77*	-0.48	-0.20	-0.20	0.17
Neuroticism	[0.94, 2.52]	[0.97, 2.55]	[-1.24, 0.29]	[-1.51, 1.11]	[-1.51, 1.12]	[-1.04, 1.43]
Openness	-0.50	-0.50	0.04	-0.02	-0.04	-0.16
Openness	[-1.29, 0.30]	[-1.28, 0.30]	[-0.76, 0.83]	[-1.31, 1.25]	[-1.33, 1.24]	[-1.40, 1.09]

Logistic Regression (SES)

Multinomial logistic regression models were built that regressed BMI category onto parental socio-economic status and adolescent individual differences. Two basic models were constructed: one that hypothesized parental SES:

\[BMI_i = b_0 + b_1(SES_i) + b_2(ID_i) + e_i\]

and an individual difference were two independent predictors of BMI, and a second that hypothesized these variables interacted with each other:

\[BMI_i = b_0 + b_1(SES_i) + b_2(ID_i) + b_3(SES_i\times ID_i) + e_i\] We iterated through all individual differences – the broad Big Five personality traits, the narrow SPI-27 traits, and cognitive ability – and tested each one independently in the model as an individual difference.

Models were estimated separately for men and women.

# end goal of wrangling is a data frame of data frames
# nested dataframes correspond to a single personality trait
# score refers to a participant's score on that trait
# we also standardize each of our variables within gender


sapa_male_trait = sapa_male %>%
  dplyr::select(-starts_with("p1"), -starts_with("p2"), -starts_with("edu")) %>%
  mutate(BMI_c = factor(BMI_c, levels = c("Normal Weight", "Underweight", "Overweight", "Obese"))) %>%
  mutate(set = ifelse(row_number() %in% train_male[,1], "train", "test")) %>%
  gather("trait_name", "trait_score", -ses, -BMI_c, -BMI, -BMI_p, -set) %>%
  group_by(trait_name, set) %>%
  mutate(trait_score = scale(trait_score)) %>%
  ungroup() %>%
  group_by(trait_name) %>%
  nest()

sapa_female_trait = sapa_female %>%
  dplyr::select(-starts_with("p1"), -starts_with("p2"), -starts_with("edu")) %>%
  mutate(BMI_c = factor(BMI_c, levels = c("Normal Weight", "Underweight", "Overweight", "Obese"))) %>%
  mutate(set = ifelse(row_number() %in% train_male[,1], "train", "test")) %>%
  gather("trait_name", "trait_score", -ses, -BMI_c, -BMI, -BMI_p, -set) %>%
  group_by(trait_name, set) %>%
  mutate(trait_score = scale(trait_score)) %>%
  ungroup() %>%
  group_by(trait_name) %>%
  nest()

Controlling for personality

To estimate the effect of socioeconomic status on BMI category, we graph the estimates of the SES slope coefficient across all logistic regression models controlling for individual differences. This presents not only the average estimate across all models (solid line), but the range of estimates – a wide range suggests that the effect of SES on BMI is sensitive to the inclusion of different individual difference measures, while a narrow range suggests that the effect of SES on BMI is persistent through personality and cognition.

Female

Parental socioeconomic status positively predicted greater likelihood of all non-normal categories (Underweight, Overweight, and Obese) compared to Normal among girls. Adolescent girls living in higher SES households were, on average, 103% less likely to be Underweight, 133% less likely, and 160% less likely to be Obese compared to low SES counterparts.

Male

Parental socioeconomic status positively predicted greater likelihood of all non-normal categories (Underweight, Overweight, and Obese) compared to Normal among boys. Adolescent boys living in higher SES households were, on average, 100% less likely to be Underweight, 126% less likely, and 144% less likely to be Obese compared to low SES counterparts.

Interaction with personality

To estimate the joint effect of socioeconomic status and individual differences on BMI category, we graph the estimates of the interaction terms of SES by individual differences by BMI category. Like before, we present the average effect (solid black line) and the 95% confidence intervals for each model.

Female

Parental socioeconomic status positively predicted greater likelihood of all non-normal categories (Underweight, Overweight, and Obese) compared to Normal among girls. Adolescent girls living in higher SES households were, on average, 101% less likely to be Underweight, 97% less likely, and 103% less likely to be Obese compared to low SES counterparts.

Male

Parental socioeconomic status positively predicted greater likelihood of all non-normal categories (Underweight, Overweight, and Obese) compared to Normal among boys. Adolescent boys living in higher SES households were, on average, -98% more likely to be Underweight, -104% more likely, and -100% more likely to be Obese compared to low SES counterparts.

Code (Female)

female_ses_only = train(BMI_c ~ ses, data = sapa_female, 
                      subset = train_female, 
                      method = "multinom",
                      na.action = "na.exclude", 
                      trControl = ctrl)


female_log = sapa_female_trait %>%
  # train models on training subset; use mulintomial logistic regression; use specific formula
  mutate(
    cov = map(data, ~train(BMI_c ~ trait_score + ses, data = ., 
                           subset = train_female, 
                           method = "multinom",
                           na.action = "na.exclude", 
                           trControl = ctrl)),
    int = map(data, ~train(BMI_c ~ trait_score*ses, data = ., 
                           subset = train_female, 
                           method = "multinom",
                           na.action = "na.exclude", 
                           trControl = ctrl))) %>%
  gather("model", "output", cov, int) %>%
  # create test data from all rows not used in training
  mutate(test_data = map(data, .f = function(x) x[-train_female, ]),
         #extract reference (true) BMI categories from test data
         test_reference = map(test_data, "BMI_c"),
         # predict categories from model output; na.pass puts NAs in any row with missing data
         predicted = map2(output, test_data, predict, na.action = "na.pass"),
         # calculate accuracy, sensitivity, specificity, etc
         confusion = map2(predicted, test_reference, confusionMatrix),
         # extract final model coefficients
         final_mod = map(output, "finalModel"),
         # tidy output for printing
         coef = map(final_mod, broom::tidy, conf.int = TRUE))

Code (Male)

male_ses_only = train(BMI_c ~ ses, data = sapa_male, 
                 subset = train_male, 
                 method = "multinom",
                 maxit= 1000,
                 na.action = "na.exclude", 
                 trControl = ctrl)

accuracy = predict(male_ses_only, type="raw", newdata=sapa_male[-train_male, ])
postResample(sapa_male[-train_male, "BMI_c"], accuracy)

male_log = sapa_male_trait %>%
  # train models on training subset; use mulintomial logistic regression; use specific formula
  mutate(
    cov = map(data, ~train(BMI_c ~ trait_score + ses, data = ., 
                               subset = train_male, 
                               method = "multinom",
                               na.action = "na.exclude", 
                               trControl = ctrl)),
    int = map(data, ~train(BMI_c ~ trait_score*ses, data = ., 
                               subset = train_male, 
                               method = "multinom",
                               na.action = "na.exclude", 
                               trControl = ctrl))) %>%
  gather("model", "output", cov, int) %>%
  # create test data from all rows not used in training
  mutate(test_data = map(data, .f = function(x) x[-train_male, ]),
         #extract reference (true) BMI categories from test data
         test_reference = map(test_data, "BMI_c"),
         # predict categories from model output; na.pass puts NAs in any row with missing data
         predicted = map2(output, test_data, predict, na.action = "na.pass"),
         # calculate accuracy, sensitivity, specificity, etc
         confusion = map2(predicted, test_reference, confusionMatrix),
         # extract final model coefficients
         final_mod = map(output, "finalModel"),
         # tidy output for printing
         coef = map(final_mod, broom::tidy, conf.int = TRUE))

Logistic Regression (Individual differences)

Female

Trait	Obese	Overweight	Underweight
Cognitive Ability	0.82	0.76*	0.98
NA	[0.61, 1.12]	[0.60, 0.97]	[0.84, 1.14]
SPI: 27 Factors
Compassion	0.94	0.99	1.17*
NA	[0.72, 1.23]	[0.78, 1.25]	[1.00, 1.37]
Irritability	1.07	1.12	1.13
NA	[0.81, 1.43]	[0.89, 1.43]	[0.97, 1.32]
Sociability	0.95	0.77*	0.67*
NA	[0.70, 1.28]	[0.61, 0.97]	[0.58, 0.78]
Well Being	0.67*	0.64*	0.66*
NA	[0.50, 0.91]	[0.49, 0.84]	[0.56, 0.77]
Sensation Seeking	0.95	1.18	1.00
NA	[0.72, 1.24]	[0.94, 1.49]	[0.87, 1.16]
Anxiety	1.07	1.04	1.13
NA	[0.79, 1.44]	[0.81, 1.34]	[0.97, 1.32]
Honesty	1.18	0.91	1.02
NA	[0.81, 1.73]	[0.70, 1.20]	[0.86, 1.20]
Industry	0.67*	0.91	0.77*
NA	[0.50, 0.90]	[0.72, 1.15]	[0.67, 0.89]
Intellect	0.97	1.12	0.94
NA	[0.75, 1.26]	[0.87, 1.43]	[0.81, 1.08]
Creativity	0.86	0.90	1.12
NA	[0.66, 1.11]	[0.71, 1.13]	[0.96, 1.30]
Impulsivity	1.00	1.16	1.20*
NA	[0.76, 1.31]	[0.91, 1.49]	[1.04, 1.38]
Attention Seeking	1.09	0.91	0.81*
NA	[0.82, 1.44]	[0.70, 1.17]	[0.70, 0.94]
Order	0.73*	0.88	0.87
NA	[0.55, 0.97]	[0.68, 1.13]	[0.75, 1.02]
Authoritarianism	1.15	1.11	0.89
NA	[0.83, 1.60]	[0.87, 1.41]	[0.77, 1.03]
Charisma	1.16	1.01	0.83*
NA	[0.89, 1.52]	[0.77, 1.33]	[0.72, 0.96]
Trust	0.79	0.83	0.74*
NA	[0.59, 1.06]	[0.64, 1.08]	[0.63, 0.85]
Humor	0.94	1.17	0.93
NA	[0.70, 1.26]	[0.89, 1.55]	[0.80, 1.07]
Emotional Expressiveness	1.03	0.79	0.86*
NA	[0.79, 1.35]	[0.62, 1.00]	[0.75, 0.99]
Art Appreciation	0.95	1.03	1.52*
NA	[0.71, 1.26]	[0.82, 1.31]	[1.26, 1.83]
Introspection	0.80	0.83	1.22*
NA	[0.62, 1.03]	[0.66, 1.06]	[1.03, 1.44]
Perfectionism	1.01	1.05	0.84*
NA	[0.76, 1.35]	[0.82, 1.34]	[0.73, 0.97]
Self Control	0.76	0.79	0.94
NA	[0.57, 1.02]	[0.61, 1.03]	[0.81, 1.10]
Conformity	0.80	1.14	0.82*
NA	[0.61, 1.05]	[0.89, 1.46]	[0.71, 0.96]
Adaptability	0.97	0.91	0.77*
NA	[0.73, 1.29]	[0.70, 1.17]	[0.66, 0.90]
Easy Goingness	1.14	1.23	1.12
NA	[0.84, 1.55]	[0.94, 1.61]	[0.95, 1.31]
Emotional Stability	0.84	0.80	0.75*
NA	[0.66, 1.08]	[0.63, 1.02]	[0.65, 0.87]
Conservatism	0.85	0.91	0.87
NA	[0.63, 1.14]	[0.71, 1.18]	[0.75, 1.02]
SPI: 5 Factors
Agreeableness	0.99	0.85	0.99
NA	[0.74, 1.31]	[0.68, 1.06]	[0.85, 1.14]
Conscientiousness	0.70*	0.76	0.73*
NA	[0.54, 0.91]	[0.58, 1.00]	[0.63, 0.85]
Extraversion	0.98	0.86	0.80*
NA	[0.73, 1.30]	[0.68, 1.09]	[0.70, 0.92]
Neuroticism	1.73*	1.21	1.34*
NA	[1.27, 2.34]	[0.94, 1.54]	[1.15, 1.56]
Openness	0.88	0.87	1.19*
NA	[0.66, 1.17]	[0.67, 1.13]	[1.03, 1.38]

Male

Trait	Obese	Overweight	Underweight
Cognitive Ability	0.93	0.84	0.85*
NA	[0.72, 1.20]	[0.63, 1.11]	[0.73, 0.99]
SPI: 27 Factors
Compassion	1.08	1.13	1.16
NA	[0.82, 1.41]	[0.86, 1.49]	[0.99, 1.36]
Irritability	1.23	1.12	0.93
NA	[0.92, 1.65]	[0.84, 1.48]	[0.79, 1.09]
Sociability	0.74*	0.99	0.77*
NA	[0.58, 0.94]	[0.74, 1.32]	[0.66, 0.89]
Well Being	0.83	0.83	0.71*
NA	[0.64, 1.08]	[0.64, 1.09]	[0.61, 0.82]
Sensation Seeking	0.83	0.89	0.81*
NA	[0.64, 1.07]	[0.67, 1.18]	[0.70, 0.95]
Anxiety	1.01	0.83	1.17*
NA	[0.78, 1.30]	[0.65, 1.06]	[1.00, 1.37]
Honesty	0.90	0.83	1.04
NA	[0.70, 1.17]	[0.62, 1.10]	[0.88, 1.22]
Industry	0.80	1.03	0.93
NA	[0.62, 1.03]	[0.78, 1.37]	[0.79, 1.08]
Intellect	1.05	1.28	0.95
NA	[0.80, 1.38]	[0.87, 1.87]	[0.81, 1.12]
Creativity	1.17	0.98	0.97
NA	[0.90, 1.52]	[0.75, 1.29]	[0.83, 1.12]
Impulsivity	1.04	1.25	1.06
NA	[0.77, 1.41]	[0.98, 1.60]	[0.91, 1.23]
Attention Seeking	1.08	0.89	0.77*
NA	[0.83, 1.41]	[0.69, 1.16]	[0.67, 0.89]
Order	0.82	0.69*	0.82*
NA	[0.63, 1.06]	[0.53, 0.90]	[0.70, 0.95]
Authoritarianism	1.28	1.43*	1.18*
NA	[0.95, 1.75]	[1.06, 1.92]	[1.01, 1.39]
Charisma	0.84	1.12	0.84*
NA	[0.65, 1.07]	[0.87, 1.46]	[0.73, 0.97]
Trust	1.06	1.04	1.05
NA	[0.80, 1.40]	[0.78, 1.38]	[0.89, 1.22]
Humor	1.23	1.15	0.86
NA	[0.91, 1.65]	[0.86, 1.54]	[0.74, 1.00]
Emotional Expressiveness	0.77*	1.00	0.91
NA	[0.59, 0.99]	[0.77, 1.30]	[0.78, 1.05]
Art Appreciation	0.91	0.86	1.07
NA	[0.71, 1.18]	[0.67, 1.09]	[0.91, 1.26]
Introspection	0.85	1.07	1.04
NA	[0.67, 1.07]	[0.82, 1.39]	[0.89, 1.21]
Perfectionism	0.83	0.73*	1.07
NA	[0.64, 1.09]	[0.56, 0.94]	[0.92, 1.24]
Self Control	0.58*	0.87	1.05
NA	[0.45, 0.75]	[0.66, 1.14]	[0.89, 1.23]
Conformity	0.96	1.40*	0.99
NA	[0.75, 1.22]	[1.04, 1.88]	[0.85, 1.15]
Adaptability	0.89	1.01	0.99
NA	[0.68, 1.16]	[0.77, 1.33]	[0.86, 1.15]
Easy Goingness	1.58*	1.51*	1.34*
NA	[1.20, 2.09]	[1.12, 2.04]	[1.15, 1.57]
Emotional Stability	0.92	0.79	0.74*
NA	[0.71, 1.20]	[0.61, 1.00]	[0.64, 0.87]
Conservatism	0.99	0.98	0.84*
NA	[0.73, 1.35]	[0.75, 1.28]	[0.73, 0.98]
SPI: 5 Factors
Agreeableness	1.17	0.88	1.13
NA	[0.89, 1.54]	[0.67, 1.16]	[0.97, 1.32]
Conscientiousness	0.75*	0.71*	0.93
NA	[0.57, 0.99]	[0.52, 0.97]	[0.80, 1.09]
Extraversion	0.80	1.03	0.76*
NA	[0.62, 1.04]	[0.78, 1.37]	[0.65, 0.89]
Neuroticism	1.18	0.86	1.11
NA	[0.91, 1.53]	[0.67, 1.10]	[0.95, 1.29]
Openness	1.04	0.81	1.00
NA	[0.80, 1.37]	[0.62, 1.07]	[0.86, 1.17]

Model accuracy

packages = c("glmnet", "caret", "ggpubr", "knitr", "kableExtra")
lapply(packages, library, character.only = TRUE)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## Loaded glmnet 4.1-1

rm(packages)

load("data/cleaned.Rdata")

We use the splits generated in the Cleaning data section of the notebook to create our training and test data. (Note that to this point, only the training data have been used in the regression models).

female_train = sapa_female[train_female, ] %>% select(BMI_p, ses, cog, contains("SPI")) %>% filter(complete.cases(.))
female_test = sapa_female[-train_female, ] %>% select(BMI_p, ses, cog, contains("SPI")) %>% filter(complete.cases(.))
female_bmi = female_train$BMI_p

male_train = sapa_male[train_male, ] %>% select(BMI_p, ses, cog, contains("SPI")) %>% filter(complete.cases(.))
male_test = sapa_male[-train_male, ] %>% select(BMI_p, ses, cog, contains("SPI")) %>% filter(complete.cases(.))
male_bmi = male_train$BMI_p

We build a function to fit the lasso models.

fit_model = function(data, outcome){
  cv_value = model.matrix(BMI_p ~ ., 
                 data = data) %>%
    cv.glmnet(x = .,
              y = outcome, 
              alpha = 1)
  model =  model.matrix(BMI_p ~ ., data = data) %>%
    glmnet(y = outcome, 
           alpha = 1, 
           lambda = cv_value$lambda.min)
  return(model)
}

We also a build a model to get the predictions in the test set from the model best fit in the training data.

pred_model = function(model, test.data){
  
  if(length(model$coefficients) == 2){
    x = test.data
    }else{ 
      x = test.data[, c("BMI_p", rownames(model$beta)[-1])]
      x <- model.matrix(BMI_p~., x)
      }
  predictions = model %>% predict(x) %>% as.vector()
  # Model performance metrics
  fit = data.frame(
    RMSE = RMSE(predictions, test.data$BMI_p),
    Rsquare = R2(predictions, test.data$BMI_p)
  )
  return(fit)
}

We fit these models separately for adolescent boys and adolescent girls

set.seed(060821)

mod1_f = lm(BMI_p ~ ses, data = female_train)
mod2_f = female_train %>%
  select(BMI_p, ses, cog) %>%
  fit_model(data = ., outcome = female_bmi)
mod3_f = female_train %>%
  select(BMI_p, ses, contains("SPI")) %>%
  select(1:7) %>%
  fit_model(data = ., outcome = female_bmi)
mod4_f = female_train %>%
  select(BMI_p, ses, contains("SPI")) %>%
  select(1:2,8:34) %>%
  fit_model(data = ., outcome = female_bmi)
mod5_f = female_train %>%
  select(BMI_p, ses, cog, contains("SPI")) %>%
  select(1:8) %>%
  fit_model(data = ., outcome = female_bmi)
mod6_f = female_train %>%
  select(BMI_p, ses, cog, contains("SPI")) %>%
  select(1:3,9:35) %>%
  fit_model(data = ., outcome = female_bmi)

female_fits = data.frame(
  vars = c(
    "SES only",
    "SES + Cog",
    "SES + Big Five",
    "SES + Narrow 27",
    "SES + Cog + Big Five",
    "SES + Cog + Narrow 27"))
female_fits$model = list(mod1_f, mod2_f, mod3_f, mod4_f, mod5_f, mod6_f)
female_fits = mutate(female_fits, fits = map(model, pred_model, test.data = female_test))

female_fits = female_fits %>%
  select(-model) %>%
  unnest(cols = c(fits)) %>%
  mutate(gender = "Adolescent Girls")

set.seed(060821)

mod1_m = lm(BMI_p ~ ses, data = male_train)
mod2_m = male_train %>%
  select(BMI_p, ses, cog) %>%
  fit_model(data = ., outcome = male_bmi)
mod3_m = male_train %>%
  select(BMI_p, ses, contains("SPI")) %>%
  select(1:7) %>%
  fit_model(data = ., outcome = male_bmi)
mod4_m = male_train %>%
  select(BMI_p, ses, contains("SPI")) %>%
  select(1:2,8:34) %>%
  fit_model(data = ., outcome = male_bmi)
mod5_m = male_train %>%
  select(BMI_p, ses, cog, contains("SPI")) %>%
  select(1:8) %>%
  fit_model(data = ., outcome = male_bmi)
mod6_m = male_train %>%
  select(BMI_p, ses, cog, contains("SPI")) %>%
  select(1:3,9:35) %>%
  fit_model(data = ., outcome = male_bmi)

male_fits = data.frame(
  vars = c(
    "SES only",
    "SES + Cog",
    "SES + Big Five",
    "SES + Narrow 27",
    "SES + Cog + Big Five",
    "SES + Cog + Narrow 27"))
male_fits$model = list(mod1_m, mod2_m, mod3_m, mod4_m, mod5_m, mod6_m)
male_fits = mutate(male_fits, fits = map(model, pred_model, test.data = male_test))

male_fits = male_fits %>%
  select(-model) %>%
  unnest(cols = c(fits)) %>%
  mutate(gender = "Adolescent Boys")

We extract the relevant information for a table.

female_fits %>%
  full_join(male_fits) %>%
  mutate(gender = str_remove(gender, "Adolescent ")) %>%
  gather(stat, value, starts_with("R")) %>%
  unite(stat, gender, stat) %>%
  spread(stat, value) %>%
  mutate(vars = factor(vars, levels = c("SES only",
                                        "SES + Cog",
                                        "SES + Big Five",
                                        "SES + Narrow 27",
                                        "SES + Cog + Big Five",
                                        "SES + Cog + Narrow 27"))) %>%
  arrange(vars) %>%
  kable(col.names = c("Model", rep(c("RMSE", "R-squared"), 2)),
        booktabs = T,
        digits = c(0,2,3,2,3)) %>%
  kable_styling() %>%
  add_header_above(c(" ", "Adolescent Boys" = 2, "Adolescent Girls" = 2))

## Joining, by = c("vars", "RMSE", "Rsquare", "gender")

	Adolescent Boys		Adolescent Girls
Model	RMSE	R-squared	RMSE	R-squared
SES only	30.09	0.020	27.02	0.031
SES + Cog	30.02	0.024	26.95	0.036
SES + Big Five	30.11	0.020	27.02	0.030
SES + Narrow 27	29.76	0.052	26.90	0.038
SES + Cog + Big Five	30.05	0.024	26.96	0.034
SES + Cog + Narrow 27	29.64	0.055	26.86	0.041